ESTADÍSTICA 2. PRÁCTICA 2.

Práctica de Clasificación Lineal: Regresión Logística

Autor: Jesús Octavio Raboso.

Los enunciados del ejercicio aparecerán en color azul. El resto de anotaciones (explicativas y personales) no seguirán ese estilo.

***Daniel Garcia Diaz (garciad@ifca.unican.es)***

***En esta práctica utilizaremos el dataset MNIST, compuesto por 60000 imágenes de train y 10000 imagenes de test (de 28x28 pixeles) correspondientes a distintas versiones digitalizadas de los dígitos 0, ..., 9.***

***Descargamos el dataset que vamos a utilizar.***

***4 archivos diferentes para crear los datasets de train (x_train->train-image e y_train->train-labels) y de test (x->t10k-images e y->t10k-labels).***

***Funciones para cargar los archivos***

***Dibujamos los 6 primeros digitos para ver el dataset que tenemos***

***El objetivo es clasificar correctamente cada una de las imágenes, utilizando el valor de los 784 pixeles. Este problema se encuandra dentro de la clasificación multi-categoría (tenemos 10 posibles clases a predicir para cada dato).***

***Para simplificar estos problemas, se suelen consideran tantos problemas binarios (predecir sí o no) independientes como categorías se tengan. Es decir, un modelo para clasificar 0 (si o no), otro para 1, etc. En este ejemplo, consideramos el dígito '9' y el objetivo es predecir si una imagen es 9 o no.***

***Como el tamaño de la muestra es muy grande (60000 imagenes), seleccionamos una muestra para entrenar; por ejemplo, los 5000 primeros dígitos.***

***Aunque no es el modelo adecuado, ya que no está acotada y puede tomar valores mucho mayores que 1 o menores que 0, construimos primero un modelo de regresión lineal para estos datos.***

***Hacemos la salida del modelo binaria considerando el umbral en 0.5 (a modo de probabilidad) y calculamos la tasa de acierto con los datos de train***

***Ahora construimos un modelo de regresión logística utilizando la función 'glm' con la familia 'family = binomial(link = "logit")'***

***Las inestabilidades numéricas se producen por la alta dimensionalidad de los datos. Algunos de los coeficientes resultan redundantes y el proceso de optimización no converge.***

***(195 not defined because of singularities).***

***Una solución a este problema es reducir la dimensionalidad del conjunto de predictores. Para ello existen técnicas eficientes (como los métodos de regularización o las Componentes Principales) que se verán más adelante. En esta práctica utilizamos un entresacado de información, considerando sólo uno de cada 20 pixeles.***

PRÁCTICA 1

***Construir un modelo de clasificación para cada dígito, y obtener una ranking de los dígitos en base a su capacidad predictiva (el primero el que mejor se predice, etc.). Para evaluar la capacidad predictiva considerando el error de test obtenido al separar aleatoriamente la muestra en un conjunto de train (n=10000) y el resto de test.***

***Nota: usar la función 'sample'.***

***Nota2: Elige un entresacado óptimo para que la capacidad predictiva no se reduzca y se reduzca la dimensión el máximo posible.***

¿Cómo elegir un entresacado óptimo?, ¿qué píxeles seleccionar? En el ejemplo propuesto, se considera la secuencia de enteros entre 1 y 784 con paso 1 y se elige uno de cada 20 píxeles. Si bien, parece demasiado aleatorio elegir unos píxeles a voleo. Por ello, antes de elegir los píxeles y dado que no podemos utilizar técnicas como métodos de regularización o componentes principales, llevaremos a cabo tres estrategias. Elegiremos la que proporcione mejor resultado.

  1. Eliminación de márgenes + píxeles aleatorios.
  2. Eliminación de márgenes + media con vecinos adyacentes + píxeles más valorados.
  3. Eliminación de márgenes + suma de píxeles + píxeles más y menos valorados.

Considerar todos los posibles entresacados es inviable.

Estrategia 1: eliminación de márgenes + pixeles aleatorios.

Para separar aleatoriamente la muestra (almacenada en la variable x_train) emplearemos la función sample de R: https://www.rdocumentation.org/packages/base/versions/3.6.2/topics/sample.

En primer lugar, separamos la muestra x_train en un conjunto de entrenamiento de 10000 imágenes y un conjunto de test que contiene el resto. Dichos conjuntos son disjuntos:

Mantendremos estos índices en TODO el ejercicio.

En el ejemplo propuesto, se emplea el entresacado seq(1,784,20) dando lugar a una selección de 40 píxeles. Por ello, consideraré que el número máximo de píxeles asumibles es 40. Es claro que, cuantos más píxeles consideremos, maś información tendremos, pero hemos de encontrar un balance entre tasa de acierto y dimensión del problema.

En las estrategias que voy a considerar, es esperable que el mejor representante contenga 40 píxeles, pero consideraremos los resultados para 10, 20, 30 y 40 píxeles.

Es de esperar que, en cada imagen, el dígito esté centrado. Por ello, elegir píxeles de los bordes de la imágen, puede resultar contraproducente. Quizá pueda resultar mejor escoger píxeles de la zona central de la imagen que, a priori, pueden resultar más significativos. Si bien, nada nos garantiza que los dígitos vayan a estar centrados; son sólo suposiciones.

Para ver si esto se cumple, consideraremos todos los dígitos del 0 al 9. Para cada uno de los dígitos, creamos un vector binario y.digit de etiquetas (1 si se corresponde con el dígito considerado; 0 en caso contrario). Construimos el dataframe data.digit cuya primera columna es el vector de etiquetas y el resto de columnas se corresponde con los valores de los pixeles correspondientes al train. Dadas las imágenes del train que contienen el dígito que estamos estudiando, elegimos las 50 primeras (probar otros valores) y ploteamos qué valor contienen sus píxeles.

Trabajaremos únicamente con la muestra que tenemos de entrenamiento y se corresponde con las imágenes de los índices idx.train. Se supone que el test no lo podemos ver.

En la mayoría de los casos, podemos eliminar los píxeles de los márgenes, pues tienen valor 0. Las imágenes nos son dadas como listas de 784 elementos cuyos índices se corresponden con los enteros entre 1 y 784. Si entendemos las imágenes como matrices $28\times28$ de píxels, podemos hacer una correspondencia entre los índices de la lista y los de la matriz que se muestra a continuación. El elemento que ocupe el índice $k$ en la lista, estará en la posición cuyo valor sea $k$ en la siguiente matriz de índices:

\begin{bmatrix} 28\cdot0 + 1 & 28\cdot0 + 2 & ... & 28\cdot0 + 28 \\ 28\cdot1 + 1 & 28\cdot1 + 2 & ... & 28\cdot1 + 28 \\ \vdots & \vdots \\ 28\cdot27 + 1 & 28\cdot27 + 1 & ... & 28\cdot27 + 28 \end{bmatrix}

Cada posición es de la forma $k = 28\cdot i + j$ donde $i =0,1,...,27$, $j=1,2,...,28$:

A la vista de las gráficas anteriores, parece que podemos eliminar las primeras 4 filas (correspondientes a los píxeles $1,2,...,28\cdot4$) y las últimas y columnas (correspondientes a los píxeles $28\cdot24+1,...,28\cdot27 +28$).

También vamos a elminar las columnas laterales. Eliminaremos los píxeles de los índices $28\cdot i + j$ para $i=4,5,6,..,23$, $j=1,2,3,4,25,26,27,28$.

Con esto, estamos eliminando 384 píxeles.

Nos quedamos con los píxels "centrales", que son 400:

Una vez seleccionados los $400$ píxeles centrales, elegiremos un número aleatorio de estos como entresacado.

Creamos una función que implementa esta estrategia. La función strategy.1 recibe como input:

Los pasos que realiza son:

Creamos una función para mostrar gráficos:

En el ejemplo propuesto, se seleccionaban 1 de cada 20 pixeles, quedando una selección de 40. Ese era el entresacado. Dados los píxeles centrales, elegiremos, al azar, 10, 20, 30, 40 pixels:

Ploteamos sólo el accuracy sobre el test comparando según número de pixels:

Observamos que para 30 y 40 píxeles los resultados van intercalándose. En lugar de emplear $M$ píxeles aleatorios, ¿podemos filtrar de algún modo?

En la eliminación de los márgenes, hemos despreciado las 4 primeras columnas, las 4 últimas columnas, las 4 primeras filas, las 4 ultimas filas.

Por ello, tenemos índices que toman valores $\{28\cdot i + j\}_{i=4,5,...,23; j=5,6,...,24}$. Quizá podamos emplear las propiedades de los índices en cuando a artimética en módulo 28.

Seleccionemos, para obtener 40, los píxeles de las columnas correspondientes a $j=11,22$.

Veamos qué tal entresecado es el resultante de quitar márgenes y quedarnos con las columnas $j=11,22$

Seleccionemos, para tener 40 píxeles de dos columnas $j=13,21$ (podríamos coger otros valores):

Comparamos:

Puesto que hay múltiples combinaciones a la hora de filtrar por columnas (o filas), escojo como representante de esta estrategia 40 pixeles aleatorios dentro de los pixeles centrales.

Estrategia 2: media de adyacentes + eliminar márgenes + píxeles más altos.

Supongamos que tenemos una matriz con el valor de los píxeles para cada imagen. Para cada una de las entradas, consideraremos la propia entrada y los elementos adyacentes en diagonal, vertical y horizontal. Dada una entrada, sumaremos el valor de todos los píxeles adyacentes a esta (que existan) y el valor de la propia entrada. El valor de la entrada central, lo sustituiremos por la media de los valores considerados.

Hemos de tener en cuenta que no todas las entradas de la matriz tienen el mismo número de vecinos adyacentes. Supongamos una matriz:

\begin{bmatrix} a_{11} & a_{12} & ... & a_{1M} \\ a_{21} & a_{22} & ... & a_{2M} \\ \vdots & \vdots & ... & \vdots \\ a_{N1} & a_{N2} & ... & a_{NM} \\ \end{bmatrix}

Por ejemplo:

Una posibilidad es implementr un algoritmo que distinga por casos: si el número de la fila o columna se sale o no de los límites, entonces considerar tales vecinos. Otra opción es añadir una columna de ceros por la derecha y la izquierda; una fila de ceros por arriba y por abajo de la matriz original, de modo que resulte:

\begin{bmatrix} 0 & 0 & 0 & ... & 0 & 0 \\ 0 &a_{11} & a_{12} & ... & a_{1M} & 0 \\ 0 &a_{21} & a_{22} & ... & a_{2M} & 0 \\ 0 &\vdots & \vdots & ... & \vdots & 0 \\ 0 &a_{N1} & a_{N2} & ... & a_{NM} & 0 \\ 0 & 0 & 0 & ... & 0 & 0 \\ \end{bmatrix}

De este modo, recorreríamos sólo la matriz central que nos interesa y no haríamos distinciones en cuanto a vecinos adyacentes. Todas las entradas tienen 8 vecinos (más ella misma).

Esto iría en contra de la tendencia habitual de mejorar algoritmos en lugar de aumentar objetos almacenados, pero no supone un gran drama en lo que pretendo hacer.

Ejemplo de aplicación:

Para llevar a cabo esta estrategia, consideraremos todos los dígitos. Para cada uno de los dígitos, buscamos 50 imágenes del train que lo contengan. Inicializamos una matriz de ceros 28x28 que irá acumulando los resultados. Para cada una de las imágenes consideradas, aplicamos la media de adyacentes y sumamos la matriz resultante a la que inicializamos anteriormente. De este modo, la matriz inicial de ceros habrá acumulado la suma de las 50 matrices resultantes de aplicar la media de adyacentes. En esa matriz, seleccionamos las N entradas (píxeles) de mayor valor que no están en los márgenes.

La función strategy.2, que implementa la estrategia 2, recibe el número $N$ de píxeles que queremos. Para cada dígito:

Podríamos primero quitar márgenes y luego aplicar la media de adyacentes. Pero, posteriormente habría que mapaear los índices da la matriz resultante de eliminar márgenes en los índices de la original ¿qué función hace eso? ¿Alguna que tenga en cuenta divisores i y restos j ?

Consideramos 10, 20, 30, 40 píxeles:

Vemos que para 30 y 40 píxeles, los resultados se intercalan. ¿Con cuál nos quedamos? Veamos el error en el train:

Estamos en la misma situación. Hagamos la media de accuracy sobre test.

El representante será:

Estrategia 3: suma + eliminar márgenes + máximos/mínimos.

En la estrategia anterior, consierábamos 50 imágenes aleatorias que se correspondían con el dígito en estudio y, entendiéndolas como una matriz, hacíamos la media de cada entrada con su adyacentes y considerábamos los $N$ píxeles con más valor que se encontraban en la zona central resultante de eliminar los márgenes.

Seguiremos un punto de vista similar. Esta vez, no haremos la media sino que sumaremos directamente las matrices de las 50 imágenes consideradas. Para la matriz resultante, consideraremos los $N_{max}$ píxeles con mayor valor y los los $N_{min}$ píxeles con menos valor que se encuentren en la zona central.

Parece que los píxeles realmente significativos son aquellos que están coloreados, es decir, que contienen un valor más alto. Veamos que ocurre sobre el train para tener un mejor criterio:

Parece que la linea más adecuada es la verde, pero la azul y la roja tampoco son despreciables. Hagamos la media sobre el test:

Comparación estrategias para selección

Nos quedamos con la primera estrategia: selección de 40 píxeles aleatorios tras elminar márgenes. Hemos de tener en cuenta la aleatoriedad.

Seleccionamos el representante de cada estrategia. Ordenamos su accuracy en el test desde el mdígito que mejor clasifica hasta el peor: (Restamos 1 ya que los dítitos empiezna en 0 y no en 1 como sí empiezan los índices en R:

Ranking estrategia 1:

Ranking estrategia 2:

Ranking estrategia 3:

Observamos que el 1 es el dígito que mejor clasifica y el 9 el que peor lo hace. Es causa de la grafía de los propios dígitos.

***PRACTICA 2: Tener en cuenta la variabilidad del error de test a la hora de construir el ranking anterior. Para ello, además de calcular la tasa de acierto para cada dígito, considerar también un "intervalo de confianza" obtenido como la dispersión (desviación típica) de 10 medidas de test obtenidas con 10 muestras aleatorias distintas. ¿Existe algún dígito que pueda predecirse significativamente mejor que los demás?***

Seleccionaos la estrategia 1. Estamos entrenando 10x10 modelos:

Nos dice que está renombrando columnas pero no nos importa, sabemos cómo están ordenados: el volumen de datos es pequeño.

Repetimos el ranking:

Observamos, al igual que en las anteriores gráficas, que el 1 se clasifica bastante mejor que el resto. Quizá pueda deberse a su propia naturaleza. La mayoría de la spersonas escriben el 1 como un segmento, de modo que sólo ocupará ciertas columnas en la matriz de la imagen o, si está inclinado, cierta diagonal. Además, las diagonales siempre irán en el mismp sentido (no creo que mucha gente escriba el 1 inclinado empezando arriba a la izquiera y acabando abajo a la derecha).

En clase se ha sigerido que, dado un modelo para cada dígito (10 modelos en total), se dividea la parte del test en 10 trozos y se haga una predicción en cada uno de ellos. Lo hacemos con la estrategia 1.

Cada columna corresponde a un dígito: la columna 1 corresponde al dígito 0; la columna 2 corresponde al dígito 1...

Calculamos la media por columnas:

Calculamos la desviación típica por columnas:

Graficamos el intervalo de confianza:

Obtenemos el ranking:

***TRABAJO EXTRA: Hasta ahora hemos utilizado la tasa de acierto como medida de validación de los clasificadores. Para ello, las predicciones probabilísticas se tiene que convertir a binarias (utilizando un umbral para al probabilidad). Sin embargo, existen otras medidas más generales que consideran el carácter probabilístico de la predicción. https://es.wikipedia.org/wiki/Curva_ROC***

roc<-roc(out,as.factor(datT[,1]))

auc(roc)

plot(roc)

Emplearemos la estrategia 2, por cambiar.

Realizaremos este apartado empleando la estrategia 2.

Fijamos dígito y cambiamos pixeles:

Dibujamos las curvas ROC:

Ya vimos que el 1 era el que mejor clasificaba, por ello el área bajo la curva es tal.

Para pixeles fijos, cambiamos el dígito

Dibujamos curvas ROC:

Comparemos 1 (el que mejor clasificaba) con el 9 que el el que peor clasifica:

Interpretación de las tablas de contingencia.

Las tablas de contingencia tienen 4 entradas. Si las interpretamos como una matriz, los elementos se correspinden con:

La primera fila (POSITIVE) son aquellos elementos que han sido clasificados como positivos.

La segunda fila (NEGATIVE) son aquellos elementos que han sido clasificados como negativos.

Por otro lado, la primera columna (POSITIVE) son los que realmente son positivos.

La primera columna (NEGATIVE) son los que realmente son negativos.

Por tanto, los verdaderos positivos son aquellos elementos que realmente son positivos y que han sido clasificados positivos.

Los falsos positivos son los elementos que han sido clasificados como positivos pero que en realidad son negativos.

Los verdaderos negativos son aquellos elementos que realmente son negativos y que han sido clasificados negativos.

Los falsos negativos son los elementos que han sido clasificados como negativos pero que en realidad son positivos.

Es por ello que al calcular el accuracy nos quedamos con la diagonal, pues continen los elementos que han sido correctamente clasificados.

En nuestro caso, NEGATIVE se corresponde con 1; POSITIVE se corresponde con 0.

Interpretación de las curvas ROC y áreas bajo la curva AUC.

Una curva ROC consiste en graficar el TRUE POSITIVE RATE (conocido como Sensitivity) en función del FALSE POSITIVE RATE (conocido como 1-specificity) para diferentes puntos de corte (cut-offs). Cada punto de la curca representa un par (sensitivity, 1-specificity) correspondiente a cierto umbral de decisión. El AUC (Area Under the ROC Curve) es un amedida que nos permite discernir cómo de bien actúan los parámetros predictores en la clasificación de las clases.

En el eje cartesiano, supongamos que nos quedamos con el cuadrado de lado 1 y altura 1 entre los puntos (0,0), (1,0), (1,1), (0,1). Las curvas empiezan en (0,0) y acaban en (1,1)

Si la prueba fuera perfecta, es decir, no hubiese solapamiento, entonces la curva (forma una esquina) solo pasa por el punto (0,1) y cubre un área total de 1.

Una curva pegada a la diagonal o, lo que es lo mismo, un AUC de 0.5, no nos indica nada sobre el modelo construido. El modelo sería equivalente a tomar decisiones aleatorias.

Si tenemos que el área es 0, esto nos indica que el modelo hace justo lo contrario a lo que queremos: clasifica todo en la clase contraria.

Observamos de nuevo cómo influye la frecuencuia de las clases. Para 20 mm, la frecuencia de la clase 1 (lluvias intensas) es previsible que sea bastante pequeña. Otro motivo por el que la frecuancia de las clases sí importa.

Por ejemplo, si tuviésemos 100 observaciones y 99 fuesen 0, el mejor modelo sería aquel que siempre nos devuelve la clase 0. Fallaríamos sólo una vez.

Condicionante.

Para GLM, hemos usado un umbral de probabilidad de 0.5 para pasar a las clases binarias. Esto también condiciona los resultado. Si eligiésemos otro umbral, seríasn distintos.